MODULE pmd1_shell_pairs

   USE md1_global_types
   USE md1_global_consts
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: pmd1_get_shell_pairs

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE pmd1_get_shell_pairs(mode,basis,vecA,vecB,box,sh_pairs)

      CHARACTER(1),    INTENT(IN) :: mode
      TYPE(md1_basis), INTENT(IN) :: basis
      REAL(REALK),     INTENT(IN) :: vecA(3), vecB(3)
      REAL(REALK),     INTENT(IN) :: box(3)
      TYPE(md1_sh_pairs), POINTER :: sh_pairs(:)
      INTEGER(INTK) :: n_pairs

      NULLIFY(sh_pairs)
      IF (ASSOCIATED(sh_pairs)) THEN
         DEALLOCATE(sh_pairs)
         NULLIFY(sh_pairs)
      END IF

      ! Make list of non-vanishing shell pairs
      ! First get n_pairs
      CALL pmd1_make_shell_pairs(mode,sh_pairs,basis,vecA,vecB,box,n_pairs)  
      ALLOCATE(sh_pairs(n_pairs))
      ! Now store them
      CALL pmd1_make_shell_pairs(mode,sh_pairs,basis,vecA,vecB,box,n_pairs)

      WRITE(IO1,*) 'Number of periodic shells      =', basis%nshell
      WRITE(IO1,*) 'Number of periodic shell pairs =', SIZE(sh_pairs)

   END SUBROUTINE pmd1_get_shell_pairs

!-------------------------------------------------------------------------------

   SUBROUTINE pmd1_make_shell_pairs(mode,sh_pairs,basis,vecA,vecB,ref_cell,np)

      IMPLICIT NONE

      CHARACTER(1),       INTENT(IN)  :: mode
      TYPE(md1_sh_pairs), POINTER, INTENT(IN) :: sh_pairs(:)
      TYPE(md1_basis),    INTENT(IN)  :: basis
      REAL(REALK),        INTENT(IN)  :: vecA(3), vecB(3)
      REAL(REALK),        INTENT(IN)  :: ref_cell(3)
      INTEGER(INTK),      INTENT(OUT) :: np

!fixme
      REAL(REALK), PARAMETER :: thresh = -2.30258d+00 * 13.0d0
      TYPE(md1_sh_pairs), POINTER :: sh
      INTEGER(INTK) :: Ish,Jsh, I,J,Jhi, IJ
      REAL(REALK)   :: A(3), B(3), M(3), P(3), PM(3), RAB(3)
      REAL(REALK)   :: ExpA, ExpB, ExpP, ExpPI, ExpAR2
      REAL(REALK)   :: R2AB, ExpKAB

      np = 0
      Ishel: DO Ish = 1, basis%nshell
         A(:) = basis%xyz(:,basis%atom(Ish))
         A(:) = A(:) + vecA(:)
!FIXME FACTOR OF 2???         Jshel: DO Jsh = 1, Ish
         Jshel: DO Jsh = 1, basis%nshell
            B(:) = basis%xyz(:,basis%atom(Jsh))
            B(:) = B(:) + vecB(:)
            M(:) = Half * (A(:) + B(:))

            IF (mode == 'I') THEN
               ! Only include M(:) INSIDE the box
               IF ( ABS(M(1)) > ref_cell(1) ) CYCLE
               IF ( ABS(M(2)) > ref_cell(2) ) CYCLE
               IF ( ABS(M(3)) > ref_cell(3) ) CYCLE
            ELSE IF (mode == 'O') THEN
               ! Only include M(:) OUTSIDE the box
               IF ( ABS(M(1)) < ref_cell(1) ) CYCLE
               IF ( ABS(M(2)) < ref_cell(2) ) CYCLE
               IF ( ABS(M(3)) < ref_cell(3) ) CYCLE
            ELSE
               CALL fmm_quit('pmd1_make_shell_pairs')
            END IF

            RAB(:) = A(:) - B(:)
            R2AB = DOT_PRODUCT(RAB,RAB)
      
            IJ = 0
            DO I = basis%exp_id(Ish), basis%exp_id(Ish) + basis%nprim(Ish) - 1
               ExpA = basis%exp(I)
               ExpAR2 = ExpA * R2AB

               Jhi = basis%exp_id(Jsh) + basis%nprim(Jsh) - 1
!FIXME FACTOR OF 2???               IF (Ish == Jsh) Jhi = I
               DO J = basis%exp_id(Jsh), Jhi
                  ExpB = basis%exp(J)
                  ExpP = ExpA + ExpB
                  ExpPI = One / ExpP
                  ExpKAB = - ExpAR2 * ExpB * ExpPI

                  IF (ExpKAB >= thresh) THEN
                     IJ = IJ +1
                     P(:) = (ExpA * A(:) + ExpB * B(:)) * ExpPI
                     PM(:) = P(:) - M(:)
                  END IF

               END DO
            END DO
      
            IF (IJ <= 0) CYCLE
            np = np + 1
            IF (.NOT. ASSOCIATED(sh_pairs)) CYCLE

            IF (np > SIZE(sh_pairs)) STOP 'fixme'
!print '(A,3F15.5)', "M        = ", M(:)
!print '(A,3F15.5)', "ref_cell = ", ref_cell(:)
!print *, "Included!"
            sh => sh_pairs(np)
            sh%I = Ish
            sh%J = Jsh
            sh%A(:) = A(:)
            sh%B(:) = B(:)
            sh%xyz(:) = M(:)
            sh%Iang = basis%L(Ish)
            sh%Jang = basis%L(Jsh)
            sh%P%lo = basis%car_AO(Ish)
            sh%Q%lo = basis%car_AO(Jsh)
            sh%P%hi = basis%car_AO(Ish) + basis%LM(basis%L(Ish))%car_nlm -1
            sh%Q%hi = basis%car_AO(Jsh) + basis%LM(basis%L(Jsh))%car_nlm -1
   
         END DO Jshel
      END DO Ishel

   END SUBROUTINE pmd1_make_shell_pairs

!-------------------------------------------------------------------------------

END MODULE pmd1_shell_pairs
